Objectif : Exploration des données et modélisations
Date début de l’analyse : 17 mai 2022
Date de la dernière modification : 23 mai 2022
# Packages
packages = c("tidyverse", "glue", "kableExtra", "knitr", "plotly", "hrbrthemes", "gridExtra", "stats", "arrow", "lme4", "tidymodels", "broom.mixed", "multilevelmod", "dotwhisker", "rAmCharts")
## Installation des packages si besoin et chargement des librairies
package.check <- lapply(packages,
FUN = function(x) {
if (!require(x, character.only = TRUE)) {
install.packages(x, dependencies = TRUE)
library(x, character.only = TRUE)}})
# Import BSO complet (2013-2020)
data <- read_csv("../data/4.process/managing-variables_2nd-part/whole_bso.csv", col_types = cols(author_position = "c"), na = c("", "NA"))
# Création de nouvelles variables : couleur de journal finale et APC fiables
data <- data %>% mutate(oa_color_finale = case_when(oa_color_article_BSO == "diamond" ~ "diamond",
oa_color_openalex == "gold" ~ "gold",
oa_color_openalex == "hybrid" ~ "hybridOA",
TRUE ~ NA_character_), #qs Maurits : catégorie "other" ou tout en NA ? Sur quelle variable se baser alors ?
amount_apc_EUR_trusted = case_when(apc_source == "openAPC_estimation_publisher" | apc_source == "openAPC_estimation_publisher_year" ~ NA_real_,
TRUE ~ amount_apc_EUR))
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, bso_classification) %>%
summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup()
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, tier) %>%
summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% mutate(tier = as.character(tier))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>%
group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>%
summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>%
mutate(is_french_CA = as.character(is_french_CA_bso_wos_openalex_single_lang))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, oa_color_finale) %>%
summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, is_covered_by_couperin) %>%
summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, bso_classification) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :
Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.
# Sélection de variables pour les modèles
bso_mvp <- data %>%
select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>%
mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.logical(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA = as.logical(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.logical(is_covered_by_couperin))
# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>%
filter(!is.na(bso_classification),
!is.na(oa_color_finale),
!is.na(is_french_CA_bso_wos_openalex_single_lang),
!is.na(amount_apc_EUR_trusted),
!is.na(is_covered_by_couperin)) %>%
filter(apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0)
À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :
apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0).Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, tier) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(tier = as.character(tier))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :
Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.
# Sélection de variables pour les modèles
bso_mvp <- data %>%
select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>%
mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.logical(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA = as.logical(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.logical(is_covered_by_couperin))
# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>%
filter(!is.na(bso_classification),
!is.na(oa_color_finale),
!is.na(is_french_CA_bso_wos_openalex_single_lang),
!is.na(amount_apc_EUR_trusted),
!is.na(is_covered_by_couperin)) %>%
filter(apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0)
À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :
apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0).Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>%
group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>%
mutate(is_french_CA = as.character(is_french_CA_bso_wos_openalex_single_lang))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :
Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.
# Sélection de variables pour les modèles
bso_mvp <- data %>%
select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>%
mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.logical(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA = as.logical(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.logical(is_covered_by_couperin))
# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>%
filter(!is.na(bso_classification),
!is.na(oa_color_finale),
!is.na(is_french_CA_bso_wos_openalex_single_lang),
!is.na(amount_apc_EUR_trusted),
!is.na(is_covered_by_couperin)) %>%
filter(apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0)
À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :
apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0).Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, oa_color_finale) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :
Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.
# Sélection de variables pour les modèles
bso_mvp <- data %>%
select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>%
mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.logical(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.logical(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.logical(is_covered_by_couperin))
# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>%
filter(!is.na(bso_classification),
!is.na(oa_color_finale),
!is.na(is_french_CA_bso_wos_openalex_single_lang),
!is.na(amount_apc_EUR_trusted),
!is.na(is_covered_by_couperin)) %>%
filter(apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0)
À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :
apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0).Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, is_covered_by_couperin) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :
Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.
# Sélection de variables pour les modèles
bso_mvp <- data %>%
select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>%
mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.logical(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.logical(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.logical(is_covered_by_couperin))
# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>%
filter(!is.na(bso_classification),
!is.na(oa_color_finale),
!is.na(is_french_CA_bso_wos_openalex_single_lang),
!is.na(amount_apc_EUR_trusted),
!is.na(is_covered_by_couperin)) %>%
filter(apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0)
À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :
apc_has_been_paid == TRUE, amount_apc_EUR_trusted > 0).Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.
graph <- bso_pop %>% ggplot(aes(x = amount_apc_EUR_trusted)) +
geom_boxplot(fill = "#112446") +
coord_flip() +
labs(x = "Montant des APC fiables", y = " ", title = "Boxplot de Y") +
theme_ipsum() +
theme(axis.title.y = element_text(size=12), axis.title.x = element_text(size=12))
#boxplot <- ggplotly(graph)
graph1 <- bso_pop %>% ggplot(aes(x = amount_apc_EUR_trusted)) +
geom_histogram(bins = 30L, fill = "#112446") +
labs(x = "Montant des APC fiables", y = "Nombre d'articles", title = "Histogramme de Y") +
theme_ipsum() +
theme(axis.title.y = element_text(size=12), axis.title.x = element_text(size=12)) +
stat_function(fun = dnorm, args = list(mean = mean(bso_pop$amount_apc_EUR_trusted), sd = sd(bso_pop$amount_apc_EUR_trusted)), col = "white")
#histo <- ggplotly(graph1)
grid.arrange(graph, graph1, ncol = 2, nrow = 1)
L’étendue (écart entre les valeurs minimales et maximales) du montant des APC est importante, avec des valeurs allant de 2.27 à 9848.4 euros. Nous constatons 2 valeurs extrêmes sur le boxplot, où les montants sont respectivement de 9848.4 et 9028.43 euros, nous filtrons donc la base sur les montants inférieurs à 7000 euros (le 3e APC le plus élevé est de 6848.67 euros) et regardons à nouveau les graphiques d’évolution des APC par année, selon les différentes variables qualitatives.
# Suppression des valeurs aberrantes
bso_pop <- bso_pop %>% filter(amount_apc_EUR_trusted < 7000) %>% mutate(is_french_CA_bso_wos_openalex_single_lang = as.logical(is_french_CA_bso_wos_openalex_single_lang))
# Graphique / discipline
# table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, bso_classification) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
# plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
# Graphique / tier
# table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, tier) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
# plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
# Graphique / french CA
# table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>%
rename(is_french_CA = is_french_CA_bso_wos_openalex_single_lang)
# plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
# Graphique / oa_color_finale
# table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, oa_color_finale) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
# plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, is_covered_by_couperin) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
# Normalité de Y
mean <- mean(bso_pop$amount_apc_EUR_trusted, na.rm = TRUE)
sd <- sd(bso_pop$amount_apc_EUR_trusted, na.rm = TRUE)
ggplot(bso_pop, aes(bso_pop$amount_apc_EUR_trusted)) +
geom_histogram(aes(y=..density..), color="black", fill = "steelblue", alpha = 0.2) +
geom_density( color="red", size = 1, adjust = 5) +
stat_function(fun = dnorm, colour = "Black", size = 1, args = list(mean = mean, sd = sd)) +
xlim(c(min(bso_pop$amount_apc_EUR_trusted)-1, max(bso_pop$amount_apc_EUR_trusted)+1))+ ggtitle("Distribution du montant des APC fiables") +
xlab("Montant des APC fiables") + ylab("Densité") +
theme_classic()
# Test statistique de normalité
library(DescTools)
JarqueBeraTest(ts(bso_pop$amount_apc_EUR_trusted), robust = FALSE, method = "chisq") # Y non normal (on rejette H0)
Jarque Bera Test
data: ts(bso_pop$amount_apc_EUR_trusted) X-squared = 47076, df = 2, p-value < 2.2e-16
Le montant des APC fiables n’est pas normalement distribué.
Pour les modélisations, nous cherchons donc à expliquer et prédire le montant des APC fiables grâce aux variables suivantes :
# Train test split
train_size <- floor(0.8 * nrow(bso_pop)) # 80% of the sample size
## set the seed to make your partition reproducible
set.seed(123)
train_ind <- sample(seq_len(nrow(bso_pop)), size = train_size)
train <- bso_pop[train_ind, ]
test <- bso_pop[-train_ind, ]
Ensuite, nous séparons notre base de données en 2 échantillons : un échantillon d’apprentissage (train) et un échantillon de test (test). Le premier contient 80% des données, soit 95691 articles et le second contient les 20% derniers, soit 23923 articles. Tous les modèles seront construits sur l’échantillon d’apprentissage, puis seront validés / testés sur l’échantillon test ; s’agissant de données inconnues (les modèles n’auront pas été entraînés dessus), les appliquer à cet échantillon permettra de tester leur robustesse et leur capacité prédictive.
# Reg linéaire
lm1 = lm(amount_apc_EUR_trusted ~ year, data = train)
summary <- tidy(lm1)
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -77428.73049 | 2396.139753 | -32.31395 | 0 |
| year | 39.27705 | 1.187785 | 33.06747 | 0 |
Chaque année, le montant des APC (fiables) augmente de 39 euros en moyenne.
# Reg linéaire
lm2 = lm(amount_apc_EUR_trusted ~ year + tier + oa_color_finale + is_covered_by_couperin + is_french_CA_bso_wos_openalex_single_lang, data = train)
summary <- tidy(lm2)
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -60191.26359 | 2256.790367 | -26.67118 | 0 |
| year | 30.85978 | 1.118668 | 27.58619 | 0 |
| tier2 | -311.08861 | 9.693970 | -32.09094 | 0 |
| tier3 | -440.85370 | 9.574222 | -46.04590 | 0 |
| tier4 | -446.89060 | 10.263922 | -43.53995 | 0 |
| oa_color_finalehybridOA | 801.69500 | 8.181967 | 97.98316 | 0 |
| is_covered_by_couperinTRUE | -96.96469 | 9.262517 | -10.46850 | 0 |
| is_french_CA_bso_wos_openalex_single_langTRUE | -67.89926 | 4.874324 | -13.92998 | 0 |
La discipline est ici la source de variations du modèle, chaque discipline dispose de son ordonnée à l’origine.
# Modèle linéaire en utilisant le package lme4
model1_spec <- linear_reg() %>%
set_engine("lmer")
# une ordonnée à l'origine qui peut varier d'une discipline à l'autre
model1_fit <- model1_spec %>%
fit(amount_apc_EUR_trusted ~ (1 | bso_classification), data = train)
# vue générale du modèle
summary <- tidy(model1_fit)
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 1495.2871 | 105.4812 | 14.17586 |
| ran_pars | bso_classification | sd__(Intercept) | 332.7945 | NA | NA |
| ran_pars | Residual | sd__Observation | 810.2233 | NA | NA |
# les variations des APC par discipline
tidy(model1_fit, effects = "ran_vals") %>% # on calcule les valeurs des effets aléatoires
mutate(estimate = 1777 + estimate, # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
ymax = estimate + 2 * `std.error`) %>%
arrange(estimate) %>%
mutate(level = as_factor(level)) %>%
ggplot(aes(x = level, y = estimate)) +
geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 1777, linetype = 2) +
coord_flip()
Par rapport au premier modèle, on ajoute l’année comme variable explicative du modèle - qui impacte donc l’angle de la pente de régression (mais celle-ci est la même par cluster, c’est-à-dire par discipline).
## Tentons un 2e modèle, dans lequel on suppose que le montant des APC augmente au cours du temps
# on commence par transformer la variable année
train <- train %>%
mutate(year = year - 2013)
model2_fit <- model1_spec %>%
fit(amount_apc_EUR_trusted ~ year + (1 | bso_classification), data = train) # l'effet de l'année sur l'APC est supposé le même quel que soit la discipline
# vue générale du modèle
summary <- tidy(model2_fit) # l'APC moyen augmente de 34,6€ par an
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 1311.20846 | 103.378135 | 12.68361 |
| fixed | NA | year | 41.67462 | 1.174069 | 35.49587 |
| ran_pars | bso_classification | sd__(Intercept) | 325.72589 | NA | NA |
| ran_pars | Residual | sd__Observation | 804.94652 | NA | NA |
# les variations des APC par discipline - ici les valeurs sont calculées pour 2013 (year = 0)
tidy(model2_fit, effects = "ran_vals") %>% # on calcule les valeurs des effets aléatoires
mutate(estimate = 1565 + estimate, # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
ymax = estimate + 2 * `std.error`) %>%
arrange(estimate) %>%
mutate(level = as_factor(level)) %>%
ggplot(aes(x = level, y = estimate)) +
geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 1565, linetype = 2) +
coord_flip()
Par rapport au deuxième modèle, on autorise la pente liée à l’année, à varier selon la discipline (more random effects).
L’année est ici la source de variations du modèle, chaque année dispose de son ordonnée à l’origine.
# Modèle linéaire en utilisant le package lme4
model1_spec <- linear_reg() %>%
set_engine("lmer")
# une ordonnée à l'origine qui peut varier d'une année à l'autre
model1_fit <- model1_spec %>%
fit(amount_apc_EUR_trusted ~ 1 + (1 | year), data = train)
# vue générale du modèle
summary <- tidy(model1_fit)
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 1771.9066 | 41.84011 | 42.34947 |
| ran_pars | year | sd__(Intercept) | 118.0731 | NA | NA |
| ran_pars | Residual | sd__Observation | 817.8890 | NA | NA |
# les variations des APC par année
tidy(model1_fit, effects = "ran_vals") %>% # on calcule les valeurs des effets aléatoires
mutate(estimate = 1772 + estimate, # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
ymax = estimate + 2 * `std.error`) %>%
arrange(estimate) %>%
mutate(level = as_factor(level)) %>%
ggplot(aes(x = level, y = estimate)) +
geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 1772, linetype = 2) +
coord_flip()
Par rapport au premier modèle, on ajoute la discipline comme variable explicative du modèle - qui impacte donc l’angle de la pente de régression (mais celle-ci est la même par cluster, c’est-à-dire par année).
## Tentons un 2e modèle, dans lequel on suppose que le montant des APC dépend de la discipline
model2_fit <- model1_spec %>%
fit(amount_apc_EUR_trusted ~ bso_classification + (1 | year), data = train) # l'effet de la discipline sur l'APC est supposé le même quelle que soit l'année
# vue générale du modèle
summary <- tidy(model2_fit)
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 1849.35886 | 43.361568 | 42.649723 |
| fixed | NA | bso_classificationChemistry | -197.16272 | 15.207234 | -12.965061 |
| fixed | NA | bso_classificationComputer and information sciences | -553.73092 | 19.908469 | -27.813838 |
| fixed | NA | bso_classificationEarth, Ecology, Energy and applied biology | -232.68578 | 9.258577 | -25.131917 |
| fixed | NA | bso_classificationEngineering | -486.09551 | 27.113104 | -17.928434 |
| fixed | NA | bso_classificationHumanities | -882.32097 | 34.474890 | -25.593148 |
| fixed | NA | bso_classificationMathematics | -925.75223 | 34.972179 | -26.471106 |
| fixed | NA | bso_classificationMedical research | -13.18786 | 6.038334 | -2.184023 |
| fixed | NA | bso_classificationPhysical sciences, Astronomy | -227.31944 | 10.420008 | -21.815669 |
| fixed | NA | bso_classificationSocial sciences | -417.68790 | 34.692059 | -12.039871 |
| ran_pars | year | sd__(Intercept) | 122.01518 | NA | NA |
| ran_pars | Residual | sd__Observation | 804.11607 | NA | NA |
# les variations des APC par année - ici les valeurs sont calculées pour la discipline "Social"
tidy(model2_fit, effects = "ran_vals") %>% # on calcule les valeurs des effets aléatoires
mutate(estimate = 1852 + estimate, # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
ymax = estimate + 2 * `std.error`) %>%
arrange(estimate) %>%
mutate(level = as_factor(level)) %>%
ggplot(aes(x = level, y = estimate)) +
geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 1852, linetype = 2) +
coord_flip()
Par rapport au deuxième modèle, on autorise la pente liée à la discipline, à varier selon l’année (more random effects).
### Hypothèses
# Test d'homogénéité de la variance (H0)
bartlett.test(amount_apc_EUR_trusted ~ interaction(year, is_covered_by_couperin), data = bso_pop)
Bartlett test of homogeneity of variances
data: amount_apc_EUR_trusted by interaction(year, is_covered_by_couperin) Bartlett’s K-squared = 1959.9, df = 15, p-value < 2.2e-16
#-->use weights arguments from lme()
# Test d'autocorrélation des résidus
# Taille des sous-populations
table(bso_pop$is_covered_by_couperin) # échantillons déséquilibrés
FALSE TRUE 102667 16947
table(bso_pop$year) # de 8184 à 26054 (peu d'obs pour 2013-2014 car suppression des NA)
2013 2014 2015 2016 2017 2018 2019 2020 8184 9537 11251 12978 14839 16943 19828 26054
# 16 clusters
#--> passer en bayésien
# Filtre sur les Y is_french_CA manquants
sample <- data %>% filter(is.na(is_french_CA_bso_wos_openalex_single_lang))
Document sous licence ouverte